4_regression

Script in series: “Associations between early life mental health indicators and sleep in adulthood”

This script performs multiple imputation and regression analyses on the dataset associated with this project.

Dependencies

We start by loading the necessary dependencies.

library(data.table)
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.1     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::between()     masks data.table::between()
## ✖ dplyr::filter()      masks stats::filter()
## ✖ dplyr::first()       masks data.table::first()
## ✖ lubridate::hour()    masks data.table::hour()
## ✖ lubridate::isoweek() masks data.table::isoweek()
## ✖ dplyr::lag()         masks stats::lag()
## ✖ dplyr::last()        masks data.table::last()
## ✖ lubridate::mday()    masks data.table::mday()
## ✖ lubridate::minute()  masks data.table::minute()
## ✖ lubridate::month()   masks data.table::month()
## ✖ lubridate::quarter() masks data.table::quarter()
## ✖ lubridate::second()  masks data.table::second()
## ✖ purrr::transpose()   masks data.table::transpose()
## ✖ lubridate::wday()    masks data.table::wday()
## ✖ lubridate::week()    masks data.table::week()
## ✖ lubridate::yday()    masks data.table::yday()
## ✖ lubridate::year()    masks data.table::year()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors

library(magrittr)

## 
## Attaching package: 'magrittr'
## 
## The following object is masked from 'package:purrr':
## 
##     set_names
## 
## The following object is masked from 'package:tidyr':
## 
##     extract

library(parallel)
library(mice)

## 
## Attaching package: 'mice'
## 
## The following object is masked from 'package:stats':
## 
##     filter
## 
## The following objects are masked from 'package:base':
## 
##     cbind, rbind

library(lmtest)

## Loading required package: zoo
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(sandwich)
library(nnet)

Loading the data

We now load the cleaned BCS70 data and prepare them for multiple imputation.

data_merged <- fread(".\\Cleaned\\data_clean.csv") %>% data.frame()

# Replacing blank values in character vectors with NA
is.na(data_merged[,]) <- data_merged[,] == ""

# Linearly rescaling all numeric variables to between 0 and 1 to aid regression model convergence
data_merged %<>% mutate_if(is.numeric, function(x) (x - min(x, na.rm = TRUE))/(max(x, na.rm = TRUE) - min(x, na.rm = TRUE)))

# Listing dependent variables for Poisson and multinomial regression, respectively
pois_list <- c("dv_sr_sleep_abn", "dv_sd_sleep_abn", "dv_pal_sleep_abn", "dv_vdb_sleep_abn", "dv_winkler_sleep_abn")
multinom_list <- c("dv_sr_sleep_abn_cat", "dv_sd_sleep_abn_cat", "dv_pal_sleep_abn_cat", "dv_vdb_sleep_abn_cat", "dv_winkler_sleep_abn_cat")

Specifying imputation methods

We now convert variables into factors as appropriate, then specify the imputation method as predictive mean matching (pmm) for all variables.

# Merging sparse categories to aid convergence
data_merged %<>% mutate(
  # Binarising e216a (mother's regular job since child's birth) to job vs no job
  e216a = case_when(
    e216a %in% c("F+Pt Job", "Ft Job", "Pt Job") ~ "Job",
    e216a == "Other" | is.na(e216a) ~ NA_character_,
    TRUE ~ e216a
  ),
  
  # Combining classes I and II, and sending "other" categories to NA, in parental social class measures
  fclrg90 = if_else(
    fclrg90 %in% c("I Professional", "II Managerial"),
    "I/II Professional/Managerial",
    fclrg90
  ),
  a0014 = case_when(
    a0014 %in% c("SC 1", "SC 2") ~ "SC 1/2",
    a0014 %in% c("Other", "Unsupported") ~ NA_character_,
    TRUE ~ a0014
  ),
  dv_sc_age_5 = case_when(
    dv_sc_age_5 %in% c("I", "II") ~ "I/II",
    dv_sc_age_5 == "Student-Volun Wk" ~ NA_character_,
    TRUE ~ dv_sc_age_5
  ),
  dv_sc_age_16 = case_when(
    dv_sc_age_16 %in% c("I", "II") ~ "I/II",
    dv_sc_age_16 %in% c("Dead", "Student") ~ NA_character_,
    TRUE ~ dv_sc_age_16
  ),
  
  # Binarising school type to comprehensive vs non-comprehensive
  BDSTYPE = case_when(
    BDSTYPE %in% c("Comprehensive", "Secondary Modern/Technical", "Scottish LEA") ~ "Comprehensive",
    BDSTYPE %in% c("Grammar", "Independent Private", "Independent Special", "LEA Special") ~ "Non-comprehensive",
    BDSTYPE == "Other" | is.na(BDSTYPE) ~ NA_character_
  ),
  
  # Binarising ameni (lack of amenities) to any lack of amenities vs none
  ameni = if_else(
    ameni %in% c("1 occasion", "2 occasions"),
    "At least 1 occasion",
    ameni
  ),
  
  # Binarising crowd (household overcrowding) to up to 1 person/room vs more than 1
  crowd = if_else(
    crowd %in% c("Over 1 to 1.5", "Over 1.5 to 2", "Over 2"),
    "Over 1",
    crowd
  )
)


# Converting binary independent/auxiliary variables
for(var in c("a0230", "a0255", "d016a", "e216a", "BDSTYPE", "emosup", "vote97", "B9SCQ17", "B9SCQ19", "B9TEN", "ameni", "crowd", "divorce", "prmnh", "sepmumbcs", "dv_father_schl", "dv_med_5", "dv_med_10", "dv_incsource", "dv_alcohol", "dv_training", "dv_participation", "dv_organisations")){
  data_merged[, var] %<>% as.factor()
}

# Converting unordered categorical independent/auxiliary variables
for(var in c("a0014", "dv_sc_age_5", "dv_sc_age_16", "dv_marital")){
  data_merged[, var] %<>% as.factor()
}

# Converting ordered categorical independent/auxiliary variables
data_merged$a0014 %<>% factor(
  ordered = TRUE,
  levels = c("SC 1", "SC 2", "SC 3 NM", "SC 3 M", "SC 4", "SC 5")
)

data_merged$a0043b %<>% factor(
  ordered = TRUE,
  levels = c("Non Smoker", "Stopped Pre-Preg", "Stopped Dur-Preg", "Ctl Smokers 1 - 4", "Ctl Smokers 5 - 14", "Ctl Smokers >= 15")
)

data_merged$brfed %<>% factor(
  ordered = TRUE,
  levels = c("Breastfed for over 1 month", "Breastfed for under 1 month", "Not breastfed")
)

data_merged$resmove %<>% factor(
  ordered = TRUE,
  levels = c("None", "1-3", "4 or more")
)

data_merged$tenure %<>% factor(
  ordered = TRUE,
  levels = c("Owned at both time points", "Owned at one time point", "Rented at both time points")
)

data_merged$fclrg90 %<>% factor(
  ordered = TRUE,
  levels = c("I/II Professional/Managerial", "IIINM Skilled non-manual", "IIIM Skilled manual", "IV Partly skilled", "V Unskilled")
)

data_merged$dv_sc_age_5 %<>% factor(
  ordered = TRUE,
  levels = c("I/II", "IIINM", "IIIM", "IV", "V")
)

data_merged$dv_sc_age_16 %<>% factor(
  ordered = TRUE,
  levels = c("I/II", "III non-manual", "III manual", "IV", "V")
)

# Converting to factors, and making "normal" the reference category for, the multinomial regression variables
for(var in multinom_list){
  data_merged[, var] %<>% as.factor() %>% relevel("Normal")
}

# Making a copy of the dataset, deleting the multinomial regression variables from the original version, and deleting the Poisson regression variables from the new version, so that they can be imputed separately
## The binary and ternary versions of each variable are highly collinear and therefore imputing them together causes chaos
data_merged_multi <- data_merged
data_merged %<>% select(-contains("abn_cat"))
data_merged_multi %<>% select(-(contains("sleep_abn") & !contains("abn_cat")))

# Creating vector of multiple imputation methods and setting method to pmm for all variables
method <- make.method(data_merged)
method[] <- "pmm"
method_multi <- make.method(data_merged_multi)
method_multi[] <- "pmm"

# Setting BCSID to not be imputed and removing it from the predictor set
method["bcsid"] <- ""
method_multi["bcsid"] <- ""
pred <- make.predictorMatrix(data_merged)
pred_multi <- make.predictorMatrix(data_merged_multi)
pred[, "bcsid"] <- 0
pred_multi[, "bcsid"] <- 0

Performing multiple imputation

We can now use the mice package to create multiple imputed datasets from the data_merged data frame, which contains missing data.

# Establishing a cluster for parallel computation to save time
## NB: THIS SECTION IS CURRENTLY SET TO CREATE AS MANY NODES AS THERE ARE CPU CORES AVAILABLE. THIS MAY CAUSE PROBLEMS IF THERE ARE OTHER TASKS RUNNING SIMULTANEOUSLY
### if you need to reduce the number, substitute `detectCores()` with e.g. `detectCores() - 2` or `6` in BOTH appearances below
cl <- makeCluster(detectCores())

# Setting random seed
clusterSetRNGStream(cl, 1858)

# Making the data frame visible in the cluster environment
clusterExport(cl, c("data_merged", "data_merged_multi", "method", "method_multi", "pred", "pred_multi"))

# Loading the mice package in the cluster environment
clusterEvalQ(cl, library(mice))

## [[1]]
## [1] "mice"      "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [7] "methods"   "base"     
## 
## [[2]]
## [1] "mice"      "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [7] "methods"   "base"     
## 
## [[3]]
## [1] "mice"      "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [7] "methods"   "base"     
## 
## [[4]]
## [1] "mice"      "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [7] "methods"   "base"     
## 
## [[5]]
## [1] "mice"      "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [7] "methods"   "base"     
## 
## [[6]]
## [1] "mice"      "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [7] "methods"   "base"     
## 
## [[7]]
## [1] "mice"      "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [7] "methods"   "base"     
## 
## [[8]]
## [1] "mice"      "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [7] "methods"   "base"     
## 
## [[9]]
## [1] "mice"      "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [7] "methods"   "base"     
## 
## [[10]]
## [1] "mice"      "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [7] "methods"   "base"     
## 
## [[11]]
## [1] "mice"      "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [7] "methods"   "base"     
## 
## [[12]]
## [1] "mice"      "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [7] "methods"   "base"     
## 
## [[13]]
## [1] "mice"      "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [7] "methods"   "base"     
## 
## [[14]]
## [1] "mice"      "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [7] "methods"   "base"     
## 
## [[15]]
## [1] "mice"      "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [7] "methods"   "base"     
## 
## [[16]]
## [1] "mice"      "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [7] "methods"   "base"     
## 
## [[17]]
## [1] "mice"      "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [7] "methods"   "base"     
## 
## [[18]]
## [1] "mice"      "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [7] "methods"   "base"     
## 
## [[19]]
## [1] "mice"      "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [7] "methods"   "base"     
## 
## [[20]]
## [1] "mice"      "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [7] "methods"   "base"

# Running multiple imputation for Poisson variables dataset with 10 imputations per core used and a maximum of 15 iterations
## This code was run with 20 cores and therefore there are 200 imputations in total
imp_pars <-
  parLapply(cl = cl, X = 1:detectCores(), fun = function(no){
    mice(data_merged, vis = "monotone", method = method, predictorMatrix = pred, m = 10, maxit = 15, printFlag = FALSE)
  })

# Merging imputations into a single mids object 
imp_merged <- imp_pars[[1]]
for(i in 2:length(imp_pars)){
  imp_merged <- mice::ibind(imp_merged, imp_pars[[i]])
}
rm(imp_pars)

# Running multiple imputation for multinomial variables dataset with 10 imputations per core used and a maximum of 15 iterations
## This code was run with 20 cores and therefore there are 200 imputations in total
imp_pars_multi <-
  parLapply(cl = cl, X = 1:detectCores(), fun = function(no){
    mice(data_merged_multi, vis = "monotone", method = method_multi, predictorMatrix = pred_multi, m = 10, maxit = 15, printFlag = FALSE)
  })
stopCluster(cl)

# Merging imputations into a single mids object 
imp_merged_multi <- imp_pars_multi[[1]]
for(i in 2:length(imp_pars_multi)){
  imp_merged_multi <- mice::ibind(imp_merged_multi, imp_pars_multi[[i]])
}
rm(imp_pars_multi)

# Plotting trace lines to verify convergence
plot(imp_merged)

plot(imp_merged_multi)

Regression

We can now perform the regression analyses.

#######
# Function : pois_regr
# Input : ind_vars, a string containing a formula which describes the independent variables in the form "ind1 + ind2 + ind3 + ..." and var_of_interest, a string containing the variable of interest (e.g. "ind1")
# Method : The function fits modified Poisson regression models to both complete cases (for the age 5 and 10 exposures only) and the multiply imputed data using each of the dependent variables in turn and the specified independent variables, obtains the risk ratios, confidence intervals and p-values for the variable of interest from the model and binds these into a data frame
# Output : A data frame containing the risk ratios, confidence intervals and p-values obtained from modified Poisson regressions of the specified independent variable against the outcome variables, controlling for the other independent variables; these are derived from both complete cases and imputed data where possible, and the minimum E-value needed to explain away the association is also included
#######
pois_regr <- function(ind_vars, var_of_interest){
  suppressWarnings({
    poisson_regression <- function(dep_var, ind_vars, var_of_interest){
      
      imp_model_fit <- with(imp_merged, coeftest(glm(as.formula(paste(dep_var, "~", ind_vars)), family = poisson(link = "log")), vcov. = sandwich))
      imp_results <- summary(pool(imp_model_fit))
      
      ind_var <- var_of_interest
      RR <- exp(imp_results[imp_results$term == var_of_interest,]$estimate)
      RR.LCI <- exp(imp_results[imp_results$term == var_of_interest,]$estimate + qnorm(0.05/2)*imp_results[imp_results$term == var_of_interest,]$std.error)
      RR.UCI <- exp(imp_results[imp_results$term == var_of_interest,]$estimate - qnorm(0.05/2)*imp_results[imp_results$term == var_of_interest,]$std.error)
      p <- imp_results[imp_results$term == var_of_interest,]$p.value
      Evalue <- if(RR.LCI < 1 && RR.UCI > 1) 1 else if(RR.UCI < 1) (1/RR.UCI) + sqrt((1/RR.UCI)*((1/RR.UCI) - 1)) else if(RR.LCI > 1) RR.LCI + sqrt(RR.LCI*(RR.LCI - 1))

      results_frame <- data.frame(dep_var, ind_var, RR, RR.LCI, RR.UCI, p, Evalue)
      
      return(results_frame)
    }
    
    results <- lapply(pois_list, function(X) poisson_regression(X, ind_vars, var_of_interest)) %>% reduce(rbind)

    return(results)
  })
}

#######
# Function : multinom_regr
# Input : ind_vars, a string containing a formula which describes the independent variables in the form "ind1 + ind2 + ind3 + ..." (the dependent variable must be a factor) and var_of_interest, a string containing the variable of interest (e.g. "ind1")
# Method : The function fits multinomial logistic regression models to both complete cases (for the age 5 and 10 exposures only) and the multiply imputed data using each of the dependent variables in turn and the specified independent variables, obtains the relative risk ratios, confidence intervals and p-values for the variable of interest from the model and binds these into a data frame
# Output : A data frame containing the relative risk ratios, confidence intervals and p-values obtained from multinomial regressions of the specified independent variable against the dependent variable, controlling for the other independent variables; these are derived from both complete cases and imputed data where possible, and the minimum E-value needed to explain away the association is also included
#######
multinom_regr <- function(ind_vars, var_of_interest){
  suppressWarnings({
    multinomial_regression <- function(dep_var, ind_vars, var_of_interest){
      
      imp_model_fit <- with(imp_merged_multi, multinom(as.formula(paste(dep_var, "~", ind_vars)), maxit = 200, trace = FALSE))
      imp_results <- summary(pool(imp_model_fit))
      
      ind_var <- var_of_interest
      group <- imp_results[imp_results$term == var_of_interest,]$y.level
      RRR <- exp(imp_results[imp_results$term == var_of_interest,]$estimate)
      RRR.LCI <- exp(imp_results[imp_results$term == var_of_interest,]$estimate + qnorm(0.05/2)*imp_results[imp_results$term == var_of_interest,]$std.error)
      RRR.UCI <- exp(imp_results[imp_results$term == var_of_interest,]$estimate - qnorm(0.05/2)*imp_results[imp_results$term == var_of_interest,]$std.error)
      p <- imp_results[imp_results$term == var_of_interest,]$p.value
      Evalue <- sapply(1:2, function(i){if(RRR.LCI[i] < 1 && RRR.UCI[i] > 1) 1 else if(RRR.UCI[i] < 1) (1/RRR.UCI[i]) + sqrt((1/RRR.UCI[i])*((1/RRR.UCI[i]) - 1)) else if(RRR.LCI[i] > 1) RRR.LCI[i] + sqrt(RRR.LCI[i]*(RRR.LCI[i] - 1))})
      
      results_frame <- data.frame(dep_var, ind_var, group, RRR, RRR.LCI, RRR.UCI, p, Evalue)
      
      return(results_frame)
    }
    
    results <- lapply(multinom_list, function(X) multinomial_regression(X, ind_vars, var_of_interest)) %>% reduce(rbind)
    
    return(results)
  })
}

# Setting random seed
set.seed(7917)

# Performing regressions against abnormal sleep

# Rutter (age 5)
rutter_5 <- pois_regr("d119 + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + tenure + prmnh + crowd + ameni + brfed", "d119")

## Original model not retained as part of coeftest object. For additional model summary information (r.squared, df, etc.), consider passing `glance.coeftest()` an object where the underlying model has been saved, i.e.`lmtest::coeftest(..., save = TRUE)`.
## This message is displayed once per session.

print(rutter_5 %>% mutate_if(is.numeric, function(x) round(x, 3)))

##                dep_var ind_var    RR RR.LCI RR.UCI     p Evalue
## 1      dv_sr_sleep_abn    d119 1.696  1.171  2.457 0.005  1.617
## 2      dv_sd_sleep_abn    d119 1.533  0.939  2.503 0.089  1.000
## 3     dv_pal_sleep_abn    d119 0.881  0.677  1.145 0.343  1.000
## 4     dv_vdb_sleep_abn    d119 1.203  0.879  1.645 0.249  1.000
## 5 dv_winkler_sleep_abn    d119 1.784  1.296  2.455 0.000  1.915

rutter_5_multi <- multinom_regr("d119 + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + tenure + prmnh + crowd + ameni + brfed", "d119")

print(rutter_5_multi %>% mutate_if(is.numeric, function(x) round(x, 3)))

##                     dep_var ind_var group   RRR RRR.LCI RRR.UCI     p Evalue
## 1       dv_sr_sleep_abn_cat    d119  Long 1.798   0.526   6.146 0.349  1.000
## 2       dv_sr_sleep_abn_cat    d119 Short 1.891   1.180   3.030 0.008  1.642
## 3       dv_sd_sleep_abn_cat    d119  Long 1.543   0.708   3.363 0.276  1.000
## 4       dv_sd_sleep_abn_cat    d119 Short 1.547   0.847   2.825 0.157  1.000
## 5      dv_pal_sleep_abn_cat    d119  Long 0.854   0.568   1.282 0.446  1.000
## 6      dv_pal_sleep_abn_cat    d119 Short 0.878   0.473   1.632 0.681  1.000
## 7      dv_vdb_sleep_abn_cat    d119  Long 1.250   0.770   2.029 0.367  1.000
## 8      dv_vdb_sleep_abn_cat    d119 Short 1.363   0.565   3.289 0.491  1.000
## 9  dv_winkler_sleep_abn_cat    d119  Long 1.172   0.572   2.400 0.665  1.000
## 10 dv_winkler_sleep_abn_cat    d119 Short 2.352   1.385   3.995 0.002  2.115

# Rutter (age 10)
rutter_10 <- pois_regr("BD3MRUTT + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + bmi + fclrg90 + tenure + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "BD3MRUTT")

print(rutter_10 %>% mutate_if(is.numeric, function(x) round(x, 3)))

##                dep_var  ind_var    RR RR.LCI RR.UCI     p Evalue
## 1      dv_sr_sleep_abn BD3MRUTT 1.688  1.048  2.718 0.032  1.271
## 2      dv_sd_sleep_abn BD3MRUTT 2.115  1.194  3.748 0.011  1.675
## 3     dv_pal_sleep_abn BD3MRUTT 1.016  0.741  1.392 0.923  1.000
## 4     dv_vdb_sleep_abn BD3MRUTT 1.324  0.898  1.952 0.158  1.000
## 5 dv_winkler_sleep_abn BD3MRUTT 1.313  0.913  1.888 0.143  1.000

rutter_10_multi <- multinom_regr("BD3MRUTT + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + bmi + fclrg90 + tenure + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "BD3MRUTT")

print(rutter_10_multi %>% mutate_if(is.numeric, function(x) round(x, 3)))

##                     dep_var  ind_var group   RRR RRR.LCI RRR.UCI     p Evalue
## 1       dv_sr_sleep_abn_cat BD3MRUTT  Long 0.927   0.190   4.536 0.926  1.000
## 2       dv_sr_sleep_abn_cat BD3MRUTT Short 2.139   1.235   3.706 0.007  1.773
## 3       dv_sd_sleep_abn_cat BD3MRUTT  Long 1.853   0.709   4.841 0.208  1.000
## 4       dv_sd_sleep_abn_cat BD3MRUTT Short 2.446   1.181   5.065 0.017  1.643
## 5      dv_pal_sleep_abn_cat BD3MRUTT  Long 1.023   0.627   1.669 0.928  1.000
## 6      dv_pal_sleep_abn_cat BD3MRUTT Short 1.149   0.538   2.457 0.720  1.000
## 7      dv_vdb_sleep_abn_cat BD3MRUTT  Long 1.619   0.852   3.076 0.143  1.000
## 8      dv_vdb_sleep_abn_cat BD3MRUTT Short 1.904   0.652   5.559 0.239  1.000
## 9  dv_winkler_sleep_abn_cat BD3MRUTT  Long 1.034   0.447   2.395 0.938  1.000
## 10 dv_winkler_sleep_abn_cat BD3MRUTT Short 1.420   0.765   2.635 0.268  1.000

# Child Development Scale
cds <- pois_regr("dv_cds_10 + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + bmi + fclrg90 + tenure + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "dv_cds_10")

print(cds %>% mutate_if(is.numeric, function(x) round(x, 3)))

##                dep_var   ind_var    RR RR.LCI RR.UCI     p Evalue
## 1      dv_sr_sleep_abn dv_cds_10 3.578  2.042  6.267 0.000  3.501
## 2      dv_sd_sleep_abn dv_cds_10 3.972  1.991  7.925 0.000  3.395
## 3     dv_pal_sleep_abn dv_cds_10 0.822  0.562  1.204 0.315  1.000
## 4     dv_vdb_sleep_abn dv_cds_10 1.114  0.741  1.676 0.604  1.000
## 5 dv_winkler_sleep_abn dv_cds_10 2.304  1.461  3.635 0.000  2.281

cds_multi <- multinom_regr("dv_cds_10 + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + bmi + fclrg90 + tenure + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "dv_cds_10")

print(cds_multi %>% mutate_if(is.numeric, function(x) round(x, 3)))

##                     dep_var   ind_var group   RRR RRR.LCI RRR.UCI     p Evalue
## 1       dv_sr_sleep_abn_cat dv_cds_10  Long 2.224   0.438  11.296 0.335  1.000
## 2       dv_sr_sleep_abn_cat dv_cds_10 Short 5.071   2.628   9.784 0.000  4.696
## 3       dv_sd_sleep_abn_cat dv_cds_10  Long 2.734   0.995   7.517 0.052  1.000
## 4       dv_sd_sleep_abn_cat dv_cds_10 Short 4.568   2.077  10.049 0.000  3.572
## 5      dv_pal_sleep_abn_cat dv_cds_10  Long 0.855   0.465   1.572 0.614  1.000
## 6      dv_pal_sleep_abn_cat dv_cds_10 Short 1.055   0.401   2.775 0.914  1.000
## 7      dv_vdb_sleep_abn_cat dv_cds_10  Long 1.279   0.619   2.643 0.506  1.000
## 8      dv_vdb_sleep_abn_cat dv_cds_10 Short 1.610   0.515   5.031 0.413  1.000
## 9  dv_winkler_sleep_abn_cat dv_cds_10  Long 0.985   0.381   2.547 0.975  1.000
## 10 dv_winkler_sleep_abn_cat dv_cds_10 Short 3.727   1.730   8.029 0.001  2.854

# Malaise Inventory
malaise <- pois_regr("BD4MAL + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + dv_bas_g + dv_med_10 + dv_sc_age_16 + bmi + fclrg90 + tenure + divorce + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "BD4MAL")

print(malaise %>% mutate_if(is.numeric, function(x) round(x, 3)))

##                dep_var ind_var    RR RR.LCI RR.UCI     p Evalue
## 1      dv_sr_sleep_abn  BD4MAL 3.602  1.839  7.055 0.000  3.081
## 2      dv_sd_sleep_abn  BD4MAL 1.575  0.818  3.034 0.176  1.000
## 3     dv_pal_sleep_abn  BD4MAL 1.158  0.787  1.703 0.457  1.000
## 4     dv_vdb_sleep_abn  BD4MAL 1.178  0.761  1.825 0.462  1.000
## 5 dv_winkler_sleep_abn  BD4MAL 1.492  0.964  2.311 0.074  1.000

malaise_multi <- multinom_regr("BD4MAL + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + dv_bas_g + dv_med_10 + dv_sc_age_16 + bmi + fclrg90 + tenure + divorce + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "BD4MAL")

print(malaise_multi %>% mutate_if(is.numeric, function(x) round(x, 3)))

##                     dep_var ind_var group   RRR RRR.LCI RRR.UCI     p Evalue
## 1       dv_sr_sleep_abn_cat  BD4MAL  Long 2.219   0.203  24.200 0.514  1.000
## 2       dv_sr_sleep_abn_cat  BD4MAL Short 5.554   2.461  12.536 0.000  4.357
## 3       dv_sd_sleep_abn_cat  BD4MAL  Long 1.492   0.321   6.937 0.610  1.000
## 4       dv_sd_sleep_abn_cat  BD4MAL Short 1.622   0.611   4.307 0.332  1.000
## 5      dv_pal_sleep_abn_cat  BD4MAL  Long 1.255   0.722   2.180 0.422  1.000
## 6      dv_pal_sleep_abn_cat  BD4MAL Short 1.367   0.551   3.393 0.501  1.000
## 7      dv_vdb_sleep_abn_cat  BD4MAL  Long 1.165   0.577   2.354 0.670  1.000
## 8      dv_vdb_sleep_abn_cat  BD4MAL Short 1.114   0.325   3.820 0.864  1.000
## 9  dv_winkler_sleep_abn_cat  BD4MAL  Long 1.482   0.663   3.311 0.338  1.000
## 10 dv_winkler_sleep_abn_cat  BD4MAL Short 1.624   0.802   3.288 0.179  1.000

# Behavioural or emotional problems
beh_emo <- pois_regr("rd6m_1 + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + dv_bas_g + dv_med_10 + dv_sc_age_16 + bmi + fclrg90 + tenure + divorce + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "rd6m_1")

print(beh_emo %>% mutate_if(is.numeric, function(x) round(x, 3)))

##                dep_var ind_var    RR RR.LCI RR.UCI     p Evalue
## 1      dv_sr_sleep_abn  rd6m_1 1.612  1.113  2.335 0.012  1.467
## 2      dv_sd_sleep_abn  rd6m_1 1.518  0.933  2.471 0.094  1.000
## 3     dv_pal_sleep_abn  rd6m_1 1.148  0.846  1.560 0.376  1.000
## 4     dv_vdb_sleep_abn  rd6m_1 0.909  0.621  1.331 0.626  1.000
## 5 dv_winkler_sleep_abn  rd6m_1 1.049  0.739  1.488 0.790  1.000

beh_emo_multi <- multinom_regr("rd6m_1 + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + dv_bas_g + dv_med_10 + dv_sc_age_16 + bmi + fclrg90 + tenure + divorce + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "rd6m_1")

print(beh_emo_multi %>% mutate_if(is.numeric, function(x) round(x, 3)))

##                     dep_var ind_var group   RRR RRR.LCI RRR.UCI     p Evalue
## 1       dv_sr_sleep_abn_cat  rd6m_1  Long 1.312   0.259   6.655 0.743  1.000
## 2       dv_sr_sleep_abn_cat  rd6m_1 Short 1.852   1.076   3.188 0.027  1.363
## 3       dv_sd_sleep_abn_cat  rd6m_1  Long 1.149   0.489   2.699 0.750  1.000
## 4       dv_sd_sleep_abn_cat  rd6m_1 Short 2.005   0.984   4.085 0.057  1.000
## 5      dv_pal_sleep_abn_cat  rd6m_1  Long 1.451   0.805   2.614 0.217  1.000
## 6      dv_pal_sleep_abn_cat  rd6m_1 Short 1.917   0.828   4.439 0.130  1.000
## 7      dv_vdb_sleep_abn_cat  rd6m_1  Long 1.016   0.494   2.091 0.965  1.000
## 8      dv_vdb_sleep_abn_cat  rd6m_1 Short 1.074   0.402   2.869 0.887  1.000
## 9  dv_winkler_sleep_abn_cat  rd6m_1  Long 0.916   0.412   2.037 0.830  1.000
## 10 dv_winkler_sleep_abn_cat  rd6m_1 Short 1.080   0.576   2.026 0.810  1.000

# Internalising behaviour in childhood
int_beh <- pois_regr("intbcsz + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + dv_bas_g + dv_med_10 + dv_sc_age_16 + bmi + fclrg90 + tenure + divorce + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "intbcsz")

print(int_beh %>% mutate_if(is.numeric, function(x) round(x, 3)))

##                dep_var ind_var    RR RR.LCI RR.UCI     p Evalue
## 1      dv_sr_sleep_abn intbcsz 1.696  1.051  2.737 0.031  1.284
## 2      dv_sd_sleep_abn intbcsz 1.508  0.912  2.493 0.110  1.000
## 3     dv_pal_sleep_abn intbcsz 1.129  0.822  1.550 0.456  1.000
## 4     dv_vdb_sleep_abn intbcsz 1.212  0.864  1.699 0.266  1.000
## 5 dv_winkler_sleep_abn intbcsz 1.189  0.848  1.668 0.316  1.000

int_beh_multi <- multinom_regr("intbcsz + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + dv_bas_g + dv_med_10 + dv_sc_age_16 + bmi + fclrg90 + tenure + divorce + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "intbcsz")

print(int_beh_multi %>% mutate_if(is.numeric, function(x) round(x, 3)))

##                     dep_var ind_var group   RRR RRR.LCI RRR.UCI     p Evalue
## 1       dv_sr_sleep_abn_cat intbcsz  Long 2.350   0.416  13.278 0.334  1.000
## 2       dv_sr_sleep_abn_cat intbcsz Short 1.930   1.080   3.448 0.027  1.375
## 3       dv_sd_sleep_abn_cat intbcsz  Long 1.134   0.489   2.633 0.769  1.000
## 4       dv_sd_sleep_abn_cat intbcsz Short 1.903   0.965   3.756 0.064  1.000
## 5      dv_pal_sleep_abn_cat intbcsz  Long 1.259   0.791   2.005 0.332  1.000
## 6      dv_pal_sleep_abn_cat intbcsz Short 1.481   0.690   3.177 0.314  1.000
## 7      dv_vdb_sleep_abn_cat intbcsz  Long 1.429   0.768   2.659 0.261  1.000
## 8      dv_vdb_sleep_abn_cat intbcsz Short 1.473   0.557   3.892 0.435  1.000
## 9  dv_winkler_sleep_abn_cat intbcsz  Long 1.096   0.516   2.330 0.811  1.000
## 10 dv_winkler_sleep_abn_cat intbcsz Short 1.298   0.706   2.386 0.401  1.000

# Externalising behaviour in childhood
ext_beh <- pois_regr("extbcsz + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + dv_bas_g + dv_med_10 + dv_sc_age_16 + bmi + fclrg90 + tenure + divorce + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "extbcsz")

print(ext_beh %>% mutate_if(is.numeric, function(x) round(x, 3)))

##                dep_var ind_var    RR RR.LCI RR.UCI     p Evalue
## 1      dv_sr_sleep_abn extbcsz 2.280  1.461  3.557 0.000  2.282
## 2      dv_sd_sleep_abn extbcsz 1.669  0.948  2.937 0.077  1.000
## 3     dv_pal_sleep_abn extbcsz 1.053  0.734  1.511 0.778  1.000
## 4     dv_vdb_sleep_abn extbcsz 1.336  0.879  2.029 0.176  1.000
## 5 dv_winkler_sleep_abn extbcsz 1.779  1.248  2.536 0.002  1.805

ext_beh_multi <- multinom_regr("extbcsz + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + dv_bas_g + dv_med_10 + dv_sc_age_16 + bmi + fclrg90 + tenure + divorce + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "extbcsz")

print(ext_beh_multi %>% mutate_if(is.numeric, function(x) round(x, 3)))

##                     dep_var ind_var group   RRR RRR.LCI RRR.UCI     p Evalue
## 1       dv_sr_sleep_abn_cat extbcsz  Long 4.837   1.151  20.334 0.032  1.567
## 2       dv_sr_sleep_abn_cat extbcsz Short 2.821   1.558   5.107 0.001  2.491
## 3       dv_sd_sleep_abn_cat extbcsz  Long 1.194   0.495   2.880 0.693  1.000
## 4       dv_sd_sleep_abn_cat extbcsz Short 2.013   0.958   4.232 0.066  1.000
## 5      dv_pal_sleep_abn_cat extbcsz  Long 1.128   0.620   2.053 0.693  1.000
## 6      dv_pal_sleep_abn_cat extbcsz Short 1.352   0.534   3.422 0.526  1.000
## 7      dv_vdb_sleep_abn_cat extbcsz  Long 1.666   0.784   3.543 0.186  1.000
## 8      dv_vdb_sleep_abn_cat extbcsz Short 1.914   0.690   5.309 0.213  1.000
## 9  dv_winkler_sleep_abn_cat extbcsz  Long 1.126   0.513   2.474 0.767  1.000
## 10 dv_winkler_sleep_abn_cat extbcsz Short 2.452   1.173   5.128 0.018  1.623

Saving results

We now write the results of the regression analyses to .csv files.

fwrite(rutter_5, ".\\Results\\rutter_5.csv")
fwrite(rutter_10, ".\\Results\\rutter_10.csv")
fwrite(cds, ".\\Results\\cds.csv")
fwrite(malaise, ".\\Results\\malaise.csv")
fwrite(beh_emo, ".\\Results\\beh_emo.csv")
fwrite(int_beh, ".\\Results\\int_beh.csv")
fwrite(ext_beh, ".\\Results\\ext_beh.csv")
fwrite(rutter_5_multi, ".\\Results\\rutter_5_multi.csv")
fwrite(rutter_10_multi, ".\\Results\\rutter_10_multi.csv")
fwrite(cds_multi, ".\\Results\\cds_multi.csv")
fwrite(malaise_multi, ".\\Results\\malaise_multi.csv")
fwrite(beh_emo_multi, ".\\Results\\beh_emo_multi.csv")
fwrite(int_beh_multi, ".\\Results\\int_beh_multi.csv")
fwrite(ext_beh_multi, ".\\Results\\ext_beh_multi.csv")

Mediation analyses (post-hoc)

We suspect that the association between early-life mental health and sleep in adulthood may be largely mediated by adult mental health. To test this, we informally perform mediation path analyses. We have already established that all seven of the exposures are associated with at least some of the adult sleep variables, but we still need to determine (a) whether childhood mental health is associated with adult mental health (here we are using Malaise and WEMWBS scores at age 42) and (b) whether controlling for these adult mental health variables attenuates the associations we found among childhood mental health and adult sleep variables.

#######
# Function : med_pois_regr
# Input : ind_vars, a string containing a formula which describes the independent variables in the form "ind1 + ind2 + ind3 + ..." and var_of_interest, a string containing the variable of interest (e.g. "ind1")
# Method : The function fits modified Poisson regression models to the multiply imputed data using each of the mediating adult mental health variables in turn and the independent variables, obtains the risk ratios, confidence intervals and p-values for the variable of interest from the model and binds these into a data frame
# Output : A data frame containing the risk ratios, confidence intervals and p-values obtained from modified Poisson regressions of the mediating adult mental health variables against the independent variables, controlling for the other independent variables
#######
med_pois_regr <- function(ind_vars, var_of_interest){
  suppressWarnings({
    mediator_poisson_regression <- function(dep_var, ind_vars, var_of_interest){
      imp_model_fit <- with(imp_merged, coeftest(glm(as.formula(paste(dep_var, "~", ind_vars)), family = poisson(link = "log")), vcov. = sandwich))
      imp_results <- summary(pool(imp_model_fit))
      
      ind_var <- var_of_interest
      impRR <- exp(imp_results[imp_results$term == var_of_interest,]$estimate)
      impRR.LCI <- exp(imp_results[imp_results$term == var_of_interest,]$estimate + qnorm(0.05/2)*imp_results[imp_results$term == var_of_interest,]$std.error)
      impRR.UCI <- exp(imp_results[imp_results$term == var_of_interest,]$estimate - qnorm(0.05/2)*imp_results[imp_results$term == var_of_interest,]$std.error)
      impP <- imp_results[imp_results$term == var_of_interest,]$p.value

      results_frame <- data.frame(dep_var, ind_var, impRR, impRR.LCI, impRR.UCI, impP)
      
      return(results_frame)
    }
    
    results <- lapply(c("BD9MAL", "BD9WEMWB"), function(X) mediator_poisson_regression(X, ind_vars, var_of_interest)) %>% reduce(rbind)

    return(results)
  })
}

# Rutter (age 5)
rutter_5_med1 <- med_pois_regr("d119 + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + tenure + prmnh + crowd + ameni + brfed", "d119")
rutter_5_med2 <- pois_regr("d119 + BD9MAL + BD9WEMWB + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + tenure + prmnh + crowd + ameni + brfed", "d119")

print(rutter_5_med1 %>% mutate_if(is.numeric, function(x) round(x, 3)))

##    dep_var ind_var impRR impRR.LCI impRR.UCI impP
## 1   BD9MAL    d119 1.818     1.537     2.151    0
## 2 BD9WEMWB    d119 0.837     0.802     0.875    0

print(rutter_5_med2 %>% mutate_if(is.numeric, function(x) round(x, 3)))

##                dep_var ind_var    RR RR.LCI RR.UCI     p Evalue
## 1      dv_sr_sleep_abn    d119 1.264  0.881  1.813 0.204  1.000
## 2      dv_sd_sleep_abn    d119 1.276  0.788  2.065 0.323  1.000
## 3     dv_pal_sleep_abn    d119 0.869  0.667  1.133 0.301  1.000
## 4     dv_vdb_sleep_abn    d119 1.149  0.841  1.571 0.384  1.000
## 5 dv_winkler_sleep_abn    d119 1.680  1.226  2.303 0.001  1.753

# Rutter (age 10)
rutter_10_med1 <- med_pois_regr("BD3MRUTT + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + bmi + fclrg90 + tenure + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "BD3MRUTT")
rutter_10_med2 <- pois_regr("BD3MRUTT + BD9MAL + BD9WEMWB + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + bmi + fclrg90 + tenure + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "BD3MRUTT")

print(rutter_10_med1 %>% mutate_if(is.numeric, function(x) round(x, 3)))

##    dep_var  ind_var impRR impRR.LCI impRR.UCI impP
## 1   BD9MAL BD3MRUTT 2.317     1.886     2.847    0
## 2 BD9WEMWB BD3MRUTT 0.764     0.724     0.807    0

print(rutter_10_med2 %>% mutate_if(is.numeric, function(x) round(x, 3)))

##                dep_var  ind_var    RR RR.LCI RR.UCI     p Evalue
## 1      dv_sr_sleep_abn BD3MRUTT 1.103  0.696  1.748 0.676      1
## 2      dv_sd_sleep_abn BD3MRUTT 1.640  0.940  2.861 0.082      1
## 3     dv_pal_sleep_abn BD3MRUTT 1.001  0.731  1.372 0.994      1
## 4     dv_vdb_sleep_abn BD3MRUTT 1.247  0.841  1.848 0.273      1
## 5 dv_winkler_sleep_abn BD3MRUTT 1.208  0.839  1.737 0.310      1

# Child Development Scale
cds_med1 <- med_pois_regr("dv_cds_10 + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + bmi + fclrg90 + tenure + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "dv_cds_10")
cds_med2 <- pois_regr("dv_cds_10 + BD9MAL + BD9WEMWB + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + bmi + fclrg90 + tenure + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "dv_cds_10")

print(cds_med1 %>% mutate_if(is.numeric, function(x) round(x, 3)))

##    dep_var   ind_var impRR impRR.LCI impRR.UCI impP
## 1   BD9MAL dv_cds_10 1.777     1.390     2.273    0
## 2 BD9WEMWB dv_cds_10 0.847     0.794     0.903    0

print(cds_med2 %>% mutate_if(is.numeric, function(x) round(x, 3)))

##                dep_var   ind_var    RR RR.LCI RR.UCI     p Evalue
## 1      dv_sr_sleep_abn dv_cds_10 2.680  1.551  4.631 0.000  2.475
## 2      dv_sd_sleep_abn dv_cds_10 3.355  1.718  6.552 0.000  2.829
## 3     dv_pal_sleep_abn dv_cds_10 0.813  0.554  1.191 0.289  1.000
## 4     dv_vdb_sleep_abn dv_cds_10 1.068  0.713  1.600 0.750  1.000
## 5 dv_winkler_sleep_abn dv_cds_10 2.189  1.389  3.449 0.001  2.124

# Malaise Inventory
malaise_med1 <- med_pois_regr("BD4MAL + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + dv_bas_g + dv_med_10 + dv_sc_age_16 + bmi + fclrg90 + tenure + divorce + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "BD4MAL")
malaise_med2 <- pois_regr("BD4MAL + BD9MAL + BD9WEMWB + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + dv_bas_g + dv_med_10 + dv_sc_age_16 + bmi + fclrg90 + tenure + divorce + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "BD4MAL")

print(malaise_med1 %>% mutate_if(is.numeric, function(x) round(x, 3)))

##    dep_var ind_var impRR impRR.LCI impRR.UCI impP
## 1   BD9MAL  BD4MAL 7.002     3.986    12.302    0
## 2 BD9WEMWB  BD4MAL 0.668     0.588     0.759    0

print(malaise_med2 %>% mutate_if(is.numeric, function(x) round(x, 3)))

##                dep_var ind_var    RR RR.LCI RR.UCI     p Evalue
## 1      dv_sr_sleep_abn  BD4MAL 1.677  0.890  3.162 0.111      1
## 2      dv_sd_sleep_abn  BD4MAL 0.997  0.505  1.967 0.993      1
## 3     dv_pal_sleep_abn  BD4MAL 1.130  0.747  1.708 0.563      1
## 4     dv_vdb_sleep_abn  BD4MAL 1.030  0.652  1.628 0.898      1
## 5 dv_winkler_sleep_abn  BD4MAL 1.270  0.799  2.019 0.314      1

# Behavioural or emotional problems
beh_emo_med1 <- med_pois_regr("rd6m_1 + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + dv_bas_g + dv_med_10 + dv_sc_age_16 + bmi + fclrg90 + tenure + divorce + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "rd6m_1")
beh_emo_med2 <- pois_regr("rd6m_1 + BD9MAL + BD9WEMWB + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + dv_bas_g + dv_med_10 + dv_sc_age_16 + bmi + fclrg90 + tenure + divorce + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "rd6m_1")

print(beh_emo_med1 %>% mutate_if(is.numeric, function(x) round(x, 3)))

##    dep_var ind_var impRR impRR.LCI impRR.UCI  impP
## 1   BD9MAL  rd6m_1 1.404     1.159     1.701 0.001
## 2 BD9WEMWB  rd6m_1 0.906     0.848     0.968 0.004

print(beh_emo_med2 %>% mutate_if(is.numeric, function(x) round(x, 3)))

##                dep_var ind_var    RR RR.LCI RR.UCI     p Evalue
## 1      dv_sr_sleep_abn  rd6m_1 1.350  0.950  1.918 0.095      1
## 2      dv_sd_sleep_abn  rd6m_1 1.372  0.846  2.226 0.201      1
## 3     dv_pal_sleep_abn  rd6m_1 1.142  0.840  1.551 0.398      1
## 4     dv_vdb_sleep_abn  rd6m_1 0.884  0.604  1.294 0.527      1
## 5 dv_winkler_sleep_abn  rd6m_1 1.011  0.713  1.433 0.952      1

# Internalising behaviour in childhood
int_beh_med1 <- med_pois_regr("intbcsz + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + dv_bas_g + dv_med_10 + dv_sc_age_16 + bmi + fclrg90 + tenure + divorce + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "intbcsz")
int_beh_med2 <- pois_regr("intbcsz + BD9MAL + BD9WEMWB + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + dv_bas_g + dv_med_10 + dv_sc_age_16 + bmi + fclrg90 + tenure + divorce + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "intbcsz")

print(int_beh_med1 %>% mutate_if(is.numeric, function(x) round(x, 3)))

##    dep_var ind_var impRR impRR.LCI impRR.UCI impP
## 1   BD9MAL intbcsz 3.451     2.863     4.160    0
## 2 BD9WEMWB intbcsz 0.743     0.705     0.784    0

print(int_beh_med2 %>% mutate_if(is.numeric, function(x) round(x, 3)))

##                dep_var ind_var    RR RR.LCI RR.UCI     p Evalue
## 1      dv_sr_sleep_abn intbcsz 0.978  0.612  1.564 0.927      1
## 2      dv_sd_sleep_abn intbcsz 1.111  0.667  1.848 0.687      1
## 3     dv_pal_sleep_abn intbcsz 1.109  0.804  1.530 0.527      1
## 4     dv_vdb_sleep_abn intbcsz 1.116  0.785  1.586 0.540      1
## 5 dv_winkler_sleep_abn intbcsz 1.061  0.753  1.493 0.736      1

# Externalising behaviour in childhood
ext_beh_med1 <- med_pois_regr("extbcsz + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + dv_bas_g + dv_med_10 + dv_sc_age_16 + bmi + fclrg90 + tenure + divorce + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "extbcsz")
ext_beh_med2 <- pois_regr("extbcsz + BD9MAL + BD9WEMWB + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + dv_bas_g + dv_med_10 + dv_sc_age_16 + bmi + fclrg90 + tenure + divorce + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "extbcsz")

print(ext_beh_med1 %>% mutate_if(is.numeric, function(x) round(x, 3)))

##    dep_var ind_var impRR impRR.LCI impRR.UCI impP
## 1   BD9MAL extbcsz 2.435     2.023     2.932    0
## 2 BD9WEMWB extbcsz 0.759     0.709     0.812    0

print(ext_beh_med2 %>% mutate_if(is.numeric, function(x) round(x, 3)))

##                dep_var ind_var    RR RR.LCI RR.UCI     p Evalue
## 1      dv_sr_sleep_abn extbcsz 1.505  0.973  2.327 0.067  1.000
## 2      dv_sd_sleep_abn extbcsz 1.300  0.731  2.312 0.373  1.000
## 3     dv_pal_sleep_abn extbcsz 1.038  0.720  1.495 0.843  1.000
## 4     dv_vdb_sleep_abn extbcsz 1.255  0.824  1.909 0.291  1.000
## 5 dv_winkler_sleep_abn extbcsz 1.644  1.155  2.340 0.006  1.578

# Exporting mediation analysis results
fwrite(rutter_5_med1, ".\\Results\\Mediation\\rutter_5_med1.csv")
fwrite(rutter_5_med2, ".\\Results\\Mediation\\rutter_5_med2.csv")
fwrite(rutter_10_med1, ".\\Results\\Mediation\\rutter_10_med1.csv")
fwrite(rutter_10_med2, ".\\Results\\Mediation\\rutter_10_med2.csv")
fwrite(cds_med1, ".\\Results\\Mediation\\cds_med1.csv")
fwrite(cds_med2, ".\\Results\\Mediation\\cds_med2.csv")
fwrite(malaise_med1, ".\\Results\\Mediation\\malaise_med1.csv")
fwrite(malaise_med2, ".\\Results\\Mediation\\malaise_med2.csv")
fwrite(beh_emo_med1, ".\\Results\\Mediation\\beh_emo_med1.csv")
fwrite(beh_emo_med2, ".\\Results\\Mediation\\beh_emo_med2.csv")
fwrite(int_beh_med1, ".\\Results\\Mediation\\int_beh_med1.csv")
fwrite(int_beh_med2, ".\\Results\\Mediation\\int_beh_med2.csv")
fwrite(ext_beh_med1, ".\\Results\\Mediation\\ext_beh_med1.csv")
fwrite(ext_beh_med2, ".\\Results\\Mediation\\ext_beh_med2.csv")